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We present a way of numerical gap estimation applicable for one- dimensional infinite uniform 
quantum systems. Using the density matrix renormalization group method for a non-uniform 
Hamiltonian, which has deformed interaction strength of j'-th bond proportional to cosh Xj, the 
uniform Hamiltonian is analyzed as a limit of A — > 0. As a consequence of the deformation, an 
excited quasi-particle is weakly bounded around the center of the system, and kept away from 
the system boundary. Therefore, insensitivity of an estimated excitation gap of the deformed 
system to the boundary allows us to have the bulk excitation gap A(A), and shift in A(A) 
from A(0) is nearly linear in A when A < 1. Efficiency of this estimation is demonstrated 
through application to the 5=1 antiferromagnetic Heisenberg chain. Combining the above 
estimation and another one obtained from the technique of convergence acceleration for finite- 
size gaps estimated by numerical diagonalizations, we conclude that the Haldane gap is in 
[0.41047905, 0.41047931]. 

KEYWORDS: DMRG, renormalization group, excitation gap, antiferromagnetic Heisenberg chain, Haldane 
gap, exact-diagonalization method 



1. Introduction 

Analysis of elementary excitations has been one of the 
central concerns in the condensed matter physics. The 
ground state of an infinitely large quantum system, that 
has a finite excitation gap, is quite different from the 
gapless systems in its correlation properties. Precise esti- 
mation of the excitation gap is therefore important, par- 
ticularly in numerical analysis of correlated systems. 

Because of limitation in computational resources, it is 
difficult to handle directly an infinite size system, but 
precise numerical analysis is possible on data for finite 
size systems. The method of finite size scaling (FSS) has 
been employed for the extrapolation of the data to the 
infinite size limit. 1,2 

As an example of the gapped system, let us consider 
a spin-5 antiferromagnetic Heisenberg chain. When 5 is 
an integer, the system has a nonzero excitation energy 
A, which is known as the Haldane gap. 3-5 In estima- 
tion of this gap, finite size corrections should be sub- 
tracted properly from the numerical data. Reliability of 
such extrapolation procedure is partially dependent on 
the maximum of available system size that is handled by 
computation resources. In the case of the 5 = 1 chain, 
the maximum at present is around 24 by use of the Lanc- 
zos diagonalization, 6-8 but it becomes thousands by use 
of the density matrix renormalization group (DMRG) 
method. 9-12 

Appropriate choice of the boundary condition is an im- 
portant procedure for the precise estimation of the exci- 
tation gap. In case of a one-dimensional system with open 
boundary conditions, reflection at the system boundary 
occasionally gives a nontrivial contribution to the kinetic 
energy of the excited quasi particle, while the particle 
itinerates in the whole system. Such a reflection effect 
can be reduced by means of a fine tuning of the boundary 



condition. In case of the 5 = 1 chain, an efficient way is to 
put an additional 5 — 1/2 spin at each end of the system, 
and to reduce the value of </ end , the coupling constant be- 
tween the 5 = 1 and the 5 = 1/2 spins at the boundary, 
compared with J, the exchange interaction inside the sys- 
tem. The value A = 0.41050(2) was reported under the 
condition J cnd = 0.5088. 13 In order to obtain a precise 
reference data for the following study in this article, we 
swept the value of J ond , and obtained a slightly smaller 
value A = 0.41047944(27) when J cnd = 0.50866, 14 as 
shown in this paper. It should be noted that this kind of 
fine tuning at the system boundary is necessary for each 
system under study. For example, if the 5 = 1 chain con- 
tains uniaxial anisotropy, the most appropriate value of 
J ond is dependent on the anisotropy parameter. 

In this article we propose a way of erasing the bound- 
ary reflection effect, by weakly confining the excited 
quasiparticle around the center of the system. For this 
purpose, we introduce the so-called hyperbolic defor- 
mation to the one-dimensional quantum Hamiltonians, 
where interaction strength between neighboring sites is 
proportional to coshAj. Here, j is the lattice index run- 
ning from — oo to oo, and A is the deformation param- 
eter. 15 Although the interaction strength becomes po- 
sition dependent, the ground state preserves a uniform 
property for any positive A. For example, the expecta- 
tion value of the bond energy of the deformed Heisen- 
berg chain is almost position independent. This uniform 
property in the ground state can be explained from the 
geometrical interpretation of the hyperbolic deformation. 

The effect of non-uniformity in the deformed Hamilto- 
nian appears in the elementary excitation. As we show 
in the following study on the deformed 5=1 Heisenberg 
chain, an excited quasiparticle is weakly attracted to the 
center of the system, where the width of the bound state 
is proportional to 1/ VA- The corresponding excitation 
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Fig. 1. Imaginary-time axis and equal-time lines in the hyper- 
bolic plane drawn inside the Poincare disc. 

gap A(A) is nearly linear in A when A <C 1. It is shown 
that the extrapolation of A (A) to the limit A — > accu- 
rately gives the Haldane gap. The obtained value is com- 
pared with another value determined by the sequence 
interval squeeze (SIS) method. 8 The precise procedure 
of the SIS method developed by Nakano and Terai is 
explained in this article and another application of this 
method is given. This technique with the exact numeri- 
cal diagonalization gives bounds for the Haldane gap. We 
will conclude definitely that the upper bound of the Hal- 
dane gap is given by lim^o A(A) much precisely than 
the SIS method. In this paper, the present best estima- 
tion of the lower bound given by the SIS method is also 
given. 

In the next section we explain the geometric back- 
ground of the hyperbolic deformation. As an example, 
we consider a deformed tight-binding model, and its 
continuum limit. In §3, we show the distribution of 
the magnetic quasiparticle under the deformation, where 
the observed shallow bound state is in accordance with 
the tight-binding picture. In §4, we perform extrapo- 
lation A — > for the estimation of the excitation gap 
lim^o A(A). As an independent estimate of the Haldane 
gap, we give the present best result by the SIS method 
in §5. Conclusions are summarized in the last section. 

2. Hyperbolic Deformation 

Real- or imaginary-time evolution of a one-dimensional 
(ID) quantum system is related to a 2D classical sys- 
tem through so called the quantum-classical correspon- 
dence. 16-18 Our aim here is to generalize the correspon- 
dence for a general case where the classical system is on 
curved 2D spaces. Let us consider a hyperbolic plane, 
which is a 2D space with constant negative curvature. 
Suppose that there is a uniform classical field on the 
hyperbolic plane, where the local action is position inde- 
pendent. Then, how does the corresponding ID quantum 
Hamiltonian look like? We consider this problem for the 
case of imaginary-time evolution. 

Figure 1 shows the 2D hyperbolic space drawn inside 
the Poincare disc. All the geodesies are represented by 
arcs, which are perpendicular to the border circle, in- 
cluding straight lines that pass through the center of the 
disc. Let us regard the vertical line as the imaginary-time 
axis. Then all the geodesies that are perpendicular to 



this imaginary-time axis can be regarded as equal- time 
curves. Suppose that the horizontal line corresponds to 
the coordinate x of the quantum system, and consider a 
quantum state |^(t)) on this line. If the classical action 
in the lower half of the hyperbolic plane is uniform, and 
if there is no symmetry breaking such as dimerization, 
the state |^(r)) is also translationally invariant. This is 
because |^(r)) is given by imaginary time boost from 
r = — co, which is mediated by the uniform action. 
Let us consider an infinitesimal evolution 

|*(t + At)) = U[At] |*(t)> (2.1) 

from t to t+At, where £/[At] represents the imaginary- 
time boost operation. Though both |*(r)) and |*(r + 
At)) are translationally invariant, W[Ar] is not. The fact 
can be seen geometrically by considering the distance be- 
tween two points (x, t) and (x, t + At) on the hyperbolic 
plane, which is an increasing function of \x\. The distance 
can be represented as (cosh^rr) At, 19 where v is a con- 
stant which is a function of the scalar curvature of the 
hyperbolic plane. If it is possible to represent U\At\ in 
the exponential form 

W[At] = cxp ^— J h(x)(coshvx) Ardx^j = exp(— AtH) , 

(2-2) 

the corresponding Hamiltonian H is also position depen- 
dent. In this case, H is written by an integral of a local 
operator h(x), and the position dependence is explicitly 
written as 

H = J (cosh vx) h(x) dx . (2.3) 

This is an example of the hyperbolic deformation of quan- 
tum Hamiltonian in the continuous ID space. If h(x) 
contains derivatives with respect to x, the form of H be- 
comes a complicated one. So, let us introduce explicit 
construction of the hyperbolic deformation starting from 
a microscopic Hamiltonian. 

Consider a way of introducing the hyperbolic deforma- 
tion to the lattice systems. We introduce lattice points at 
x = aj, where a is the lattice constant and j is the lattice 
index, which runs in a finite range from —R/a to R/a. 
The constant R we have introduced satisfies R/a ^> 1, 
and specifies the size of the system, which ensures a nu- 
merical cutoff. Then we have a relation 

cosh^x = cosh vaj = cosh Aj , (2.4) 

where A = va is the deformation parameter that we have 
used. A discrete analogue of H in eq. (2.3) is then given 
by the following lattice Hamiltonian 

H i X ) = J2 cosh hj+i + J2 cosh A(j - §) g } , (2.5) 

3 3 

where ft • - +1 represents the neighboring interaction, 
the on-site ones. Another possible choice of the discrete 
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Hamiltonian is 



H(X) = ^coshAj h j% 



cosh A j + cosh A (j — 1) 



(2.6) 



= ^2 cosh X -j h j,j+i + cosh 2 cosh A(j - i) g 3 , 

j j 

where the coefficient of the on-site term is different from 
eq. (2.5). Since we chiefly investigate small A region, this 
difference is not conspicuous. 20 We therefore choose if (A) 
in the form of eq. (2.5) in the following, ft should be 
noted that the Hamiltonian H(X — 0) is translationally 
invariant. 

As an example of the ID lattice systems, let us consider 
a non- interacting tight-binding model. The deformed 
Hamiltonian is written as 



H TB (X) = —t ^ cosh Xj ( c\c 



■j+i + c j+i c j 



) 



- M ^coshA(i-i)ct Cj , (2.7) 



where t represents the hopping amplitude, and /x the 
chemical potential. The operators, cj and Cp appear- 
ing in eq. (2.7) are fermion creation and annihilation 
operators. 21 Since there is no interaction, all the eigen- 
states can be constructed from one-particle wave func- 
tions j = (j| v E'), where \^f) is a 1-particle eigenstate 
and (j\ is defined as {0\cj. The wave function ^- of the 
stationary state satisfies the Schrodingcr equation 



E tyj = -t cosh Xj - t cosh X(j - l) ^ 
- M coshA(j - §) . 



(2.8) 



The one-particle ground state energy E becomes zero 
when —fJ,/t = 2cosh(A/2), and the corresponding wave 
function W - becomes a constant of j. If fi is smaller 
than — 2t cosh(A/2), the ground state wave function is 
bounded around the origin j = 0. This kind of bound 
state is also observed for one particle excitation in many 
body problem, as we will see in the next section. 

Let us check the continuum limit of eq. (2.8). Substi- 
tuting the relations x — aj , A = av and the correspon- 
dence 



= *(aj) = 
to eq. (2.8), we obtain the relation 



(2.9) 



-tcoshv(x — ^) coshz/^ [^(x + a) + ty(x — a)] 
-tsinh^(a; — — ) sinhi/- [ty(x + a) — $>(x — a)] 



— y(iCOsh^(a; — — )^(x) 



(2.10) 



after some algebra. Expressing the hopping amplitude as 
t = h 2 / (2ma 2 ) , chemical potential as \i = —U — 2t, and 



taking the limit a 



0, we obtain a differential equation 



h2 9 d rr , 

- — — — cosh vx — — V U cosh vx 
2m ox ox 



*(z) 



(2.11) 

The first term in the parenthesis of the r.h.s. is the de- 
formed kinetic energy, and the second term is a kind of 
trapping potential when U > 0. 

If we consider the imaginary-time dependence of the 
wave function, the Lagrangian which draws eq. (2.11) 
from the stationary condition is given by 



■JL 
dt 



<]> + cosh vx 



h 2 d^* 
2m dx dx 



+ 17**1' 



(2.12) 



for 4'(x,t), where we have introduced the letter t for 
the imaginary-time variable, and where we have used 
the unit that satisfies h = 1. Note that the time- like 
variable t in eq. (2.1) is related to t by the relation 
(cosh^x)dt = dr, and in the x-t plane the Lagrangian 
can be represented as 



= cosh vx 



T » d T n 2 d-v* a* ^ ' 

or 2m Ox ox 



(2.13) 



for ^(x,t). The action of the system is given by 



S = 



J £'(**, T tf *, 8 X ^\ % 9 T *, X *) drdx 



y*—V + h(x) 

OT 



(cosh vx) drdx. (2.14) 



This action is actually obtained by identifying \I>(a;, r) as 
a field operator and deriving the path-integral formalism 
starting from eq. (2.1) with a local Hamiltonian 

h 2 d** a* 

h ^'^^ + m " > - < 2 - 15 > 

We note that a local deformation of the measure in the 
action gives the hyperbolic deformation. 

3. Excitation of the S — 1 Heisenberg Chain 

We consider the S = 1 antifcrromagnetic Heisenberg 
chain as an example of the ID many body systems. The 
system has finite magnetic excitation energy, which is 
known as the Haldane gap. 3-5 In numerical analyses to 
obtain the gap of an open-boundary system, there is a 
custom to put S = 1/2 spins at both ends of the system, 
in order to avoid the quasi degeneracy in the low-energy 
states. 13 ' 22 The Hamiltonian of the open-boundary S = 
1 chain is represented as 



JV-l 

h(o) = j Yl Wi 

j=-N+l 

+ ^end ( S L ' S_jv+1 + S N ■ S R ) 



(3.1) 



which includes M = 2N numbers of S = 1 spins from 
S-iv+i to S N , and the boundary 5=1/2 spins s L and 
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s R . Thus there are 2N + 2 spins in total. We count the 
number of S = 1 spins M as the size of the system. 
The parameter J > represents the antiferromagnetic 
exchange coupling between neighboring S = 1 spins S^- 
and S - +1 , and J cnd > is the coupling at the boundary 
between s L and S_ N+1 and also between S N and s R . 
Throughout this article we take J as the unit of the en- 
ergy, and use the parameterization J ond = J = 1 unless 
the value of J ond is specified. 

We introduce the hyperbolic deformation for this sys- 
tem. The deformed Hamiltonian is represented as 

N-l 

H(X) = J ^ coshAj • S j+1 (3.2) 

j=-N+l 

+ Jcnd coshAiV (s L • S_ N+1 + S N ■ S R ) . 

When A = 0, eq. (3.2) becomes eq. (3.1). Occasionally it 
is convenient to treat a system that contains odd number 
of spins, so that one of the 5 = 1 spin is just at the 
center of the system. In order to satisfy the condition, we 
introduce another type of the deformed system described 
by the Hamiltonian 

JV-2 

h(\) = j coshA (i-5)V s ;+i ( 3 - 3 ) 

j=-N+l 

+J cnd cosh X(N - \) (s L • S_ N+1 + S N _i ■ s R ) , 

where there are 2N + 1 spins in total. In this case the 
system size, which is the number of S = 1 spins, is M — 
27V- 1. 

3. 1 Gap estimation for the undeformed system 

We first estimate the value of the Haldane gap A for 
undeformed systems A = 0, in order to get reference data 
for the later study under deformation A > 0. The exci- 
tation energy from the ground state is calculated by the 
DMRG method, 9-12 as a function of J ond , the system 
size M = 2N (or 2N — 1), and the number of states kept 
m, which is increased up to 180. From the various val- 
ues of J ond for which we have performed calculations, we 
show the data for two typical values J ond = 0.5088 and 
J ond = 0.50866 in Table I. The former value is used in 
a literature 13 and the latter is an optimized one in this 
work. 

The data for J ond = 0.5088 in Table I are calculated 
under the same conditions as them 13 except for a read- 
justment of the energy origin. For precise determination 
of the lowest eigenvalue, we shift the origin of the en- 
ergy so that the ground-state energy becomes nearly 
zero. This energy shift is realized by the following pro- 
cess. First we obtain the ground state |^o) diagonaliz- 
ing the Hamiltonian in eq. (3.2) or eq. (3.3), and calcu- 
late the nearest neighbor correlation function w i f+1 = 
(^qISj • S i+1 |^ ). We then replace the neighboring inter- 
action Sj • S i+1 in the Hamiltonian by S i ■ S i+1 — w i 
with an identity operator /, and perform the same sub- 
traction also for the boundary terms where s L and s R are 
involved. This subtraction can be performed successively 
when one constructs the renormalized Hamiltonians H T 



Table I. Calculated excitation energy A M for each system with 
the size M, which is the number of S = 1 spins. The integer m is 
the number of states kept. The interaction strength at the system 
boundary is J cn d = 0.5088 (the upper series from M = 27V = 100 
to 160) or J en( j = 0.50866 (the lower series). Truncation errors 
1 — Pcis{ m ) an d 1 — ^Ex( m ) are a l so shown for the ground and 
the 1st excited states, respectively. (See text.) 



M 


m 


M 


1 - PGs(m) 


1 - P EX (m) 


100 


120 


4104951946 


5.196E-12 


7.724E-10 




140 


41 04949559 


1.445E-12 


3 009E-10 




160 


0.4104948683 


3.579E-13 


1.556E-10 




180 


0.4104948265 


4.170E-14 


6.870E-11 


120 


120 


0.4104929134 


5.194E-12 


8.192E-10 




140 


0.4104926043 


1.446E-12 


3.237E-10 




160 


0.4104924896 


3.769E-13 


1.657E-10 




180 


0.4104924345 


4.125E-14 


7.301E-11 


140 


120 


0.4104912749 


5.197E-12 


8.539E-10 




140 


0.4104908951 


1.445E-12 


3.300E-10 




160 


0.4104907533 


3.678E-13 


1.732E-10 




180 


0.4104906847 


4.270E-14 


7.615E-11 


160 


120 


0.4104900535 


5.193E-12 


8.770E-10 




140 


0.4104896021 


1.443E-12 


3.391E-10 




160 


0.4104894323 


3.695E-13 


1.799E-10 




180 


0.4104893502 


5.500E-14 


7.854E-11 


100 


160 


0.4104803729 


3.657E-13 


1.562E-10 




180 


0.4104803310 


3.373E-13 


6.875E-11 


120 


160 


0.4104802346 


3.692E-13 


1.663E-10 




180 


0.4104801793 


4.313E-14 


7.306E-11 


140 


160 


0.4104801400 


3.904E-13 


1.738E-10 




180 


0.4104800712 


4.351E-14 


7.620E-11 


160 


160 


0.4104800736 


3.690E-13 


1.798E-10 




180 


0.4104799915 


6.541E-14 


7.858E-11 



and H K for the left and the right block of the system dur- 
ing the finite-size sweeping process. The ground-state en- 
ergy of the shifted Hamiltonian thus obtained is nearly 
zero. The total amount of the energy shift can be ob- 
tained from w i i+1 . It should be noted that the above en- 
ergy shift process is important for the large-scale system, 
where the ground-state energy becomes a big number. In 
the same manner, we have to use the shifted Hamiltonian 
when we consider the hyperbolically deformed system 
with A > 0, where the absolute value of the ground-state 
energy increases exponentially with the system size. 

Truncation errors introduced to the ground state 
1 — Posi 171 ), which are listed in Table I, are calculated 
by the following standard procedures in DMRG calcula- 
tion. After sufficient numbers of finite size sweeping, we 
obtain the optimized variational ground state 

l*>= E hs§s*e r \ti)\sg)\s?)\t r ), (3.4) 

where and |£ r ) represent relevant block spin state for 
the left and the right blocks, respectively, that take at 
most m numbers of states. Creating the reduced density 
matrix for the left half of the system 

PttSSitisg' = J2 Ks„ z sf^^'sf ir > (3-5) 

and diagonalizing it to obtain eigen values w a , where we 
assume the ascending order for w a . The truncation error 
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Fig. 2. System size dependence of the excitation gap A M listed 
in Table I. Only data for m = 160 and m = 180 is given. 



Fig. 3. Nearest neighbor spin correlation function (S^S^ +1 ) of 
the singlet ground state when the system size is M = 2N — 1 = 
101. 



is then calculated as 

1 - Pcsfa) = 1 - ■ (3.6) 

a=l 

The truncation error for the lowest excited state 
1 — PEx{m) is also calculated in the same manner, 
using the optimized ground-state in the subspace where 
the total number of S z is equal to 1. 

Figure 2 shows the system size dependence of the ex- 
citation gap, which is given as a list in Table I. For 
both cases, J cnd = 0.5088 and J cnd = 0.50866, the 
gap is nearly linear in 1/M = 1/(2N). From the nu- 
merical data under the condition to = 160, we ob- 
tain A = 0.41047970(1) when J cnd = 0.5088 and A = 
0.41047968(2) when J cnd = 0.50866 by use of the fit- 
ting with respect to second order polynomials, where the 
numbers shown in the parenthesis are the mean-square 
fitting error. These two values of the gap seem to be 
consistent, however, the gap is sensitive to the value 
of m. From the data when we keep to = 180 states, 
we obtain A = 0.41047947(1) when J ond = 0.5088 and 
A = 0.41047944(1) when J cnd = 0.50866. So far as we 
have calculated, the estimated value of A is always a de- 
creasing function of m. This tendency is also discussed in 
the last section. Indeed, when J end = 0.5088, the values 
are A = 0.41048145(9) for m = 120, A = 0.41048019(4) 
for m = 140, A = 0.41047970(1) for m = 160, and 
A = 0.41047947(1) for to — 180. In this parameter range 
of to, change in A as a function of to is bigger than the 
error in the last digit of the above estimation. Thus the 
lowest value obtained so far can be regarded as the up- 
per bound for A. We therefore use the smallest value 
A = 0.41047944 when m = 180 and J cnd = 0.50866 as 
the better estimation of the upper bound of A than the 
estimation of to = 180 and J end = 0.5088. 

Note that the condition J ond = 0.50866 is valid only 
for the Hamiltonians shown in Eqs. (3.2) and (3.3). If we 
need the excitation gap for a variety of S = 1 spin chains, 
which contain anisotropy and biquadratic terms, we have 
to find another appropriate value of J cnd for each cases. 



3.2 Bounded excitation when A > 

Let us observe the ground state and the elementary ex- 
citation of the deformed chain, which is described by the 
Hamiltonian in eq. (3.2) or eq. (3.3) with A > 0. Figure 3 
shows the nearest-neighbor spin correlation functions be- 
tween S = 1 spins (SfSf +1 ) of the singlet ground state 
when the system size is M = 2N — 1 = 101. In Fig. 3, 
we look only at the boundary of the system, where local 
fluctuation of this correlation function is prominent. It is 
known that the hyperbolic deformation has an effect of 
decreasing the correlation length £. 15 In the S — 1 anti- 
ferromagnetic Hcisenberg spin chain, £ is of the order of 
unity already at A = 0, and thus the effect of hyperbolic 
deformation is not conspicuous in long-range correlation 
functions, as long as the ground state is concerned. As 
displayed in Fig. 3, the short-range correlation function 
is also not affected by the deformation in this parameter 
range of A. Although the interaction strength is posi- 
tion dependent, the spin correlation function is almost 
uniform inside the system. Thus we may say that the 
vacuum of the quasiparticle excitation is kept fixed in its 
internal structure against the deformation. This behav- 
ior is favorable to have a good convergence in the energy 
gap of the elementary excitation. 

Next, we observe the magnetic excitation. Figure 4 
shows the Z-componcnt of the local spin polarization 
(Sf) calculated for the first excited state, which is the 
lowest-energy state in the subspace where the total S z of 
the system is unity. The polarization (Sf) is positive ev- 
erywhere, unless one chooses an extremely large A. Thus 
it is possible to regard (Sf) as the distribution probabil- 
ity of the excited magnetic quasiparticle. The quasipar- 
ticle is bounded around the center of the system when 
A = 0.05 and 0.10, in contrast to the unbounded case 
when A = 0. 

The observed quasiparticle distribution in Fig. 4 is 
close to the Gaussian distribution around the origin 
j = 0, when A = 0.05 or 0.10. In order to quantify the 
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Fig. 4. Spin polarization (Sj) of the first excited state. 



Fig. 6. Comparison between (n .) of the tight-binding model and 
(S?) of the Heisenberg chain. 
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Fig. 7. Linear dependence between Atb and Ahc- 



distribution width, we introduce 



Ax = 



N-l 

£ 

j=-N+l 



(3.7) 



for the cases when the system size is odd. Figure 5 
shows the value of Ax calculated for the 101-site sys- 
tem. The dotted line is the fitting for those Ax in the 
range A" 1 / 2 < 4, equivalently A > 0.04. The distribution 
width Ax is proportional to I/a/A in this parameter re- 
gion of A, where Ax is almost independent of the system 
size. We have confirmed that the relation Ax oc l/VA 
holds for l/VA < 20, equivalently A > 0.0025, when the 
system size is M = 2N = 1000. If A becomes too small 
for a fixed system size, Ax deviates from the fitting line 
due to the finite size effect. The deviation of Ax from 
the 1 / VA behavior suggests breaking the confinement of 
the quasiparticle. 

To speak qualitatively, the observed A-dependence of 
Ax is consistent with the effective one-particle potential 



U(j) = J coshAj 



(3.8) 



in the neighborhood of the origin j = 0. Note that a 
large system size M ensures existence of a finite central 
region of the chain satisfying Xj <C 1 for small but finite 



A. Thus in a final simulation, we need to have an enough 
size M <~ 1000, which is tractable by the DMRG method 
at present. 

It is possible to interpret the relation Ax\/A ~ const. 
as the quantum mechanical uncertainty for the excited 
quasiparticle under this harmonic potential. The quasi- 
particle distribution in Fig. 4 can be explained quan- 
titatively by means of the tight-binding model given 
in eq. (2.7). Let us consider the occupation number 
(nj = (cjCj) for the lowest energy one-particle state. It is 
possible to obtain a good approximation of (Sj) by (n^) if 
we choose the parametrization t = 1 and [i = 2.033. Fig- 
ure 6 shows the correspondence between (Sj) and (n^), 
where A HC represents the deformation parameter for the 
spin chain and A TB that for the tight-binding model. The 
value of A TB is determined so as to have a best fit of (S?) 
and {rij). Figure 7 shows the relation between A TB and 
A HC . We can see that a simple relation A TB = 1.324A HC 
holds, where the proportional constant 1.324 gives the 
correction to the qualitative description in eq. (3.8). 

4. Gap Estimation by Extrapolation in A 

We have observed that the magnetic excitation of 
the hyperbolically deformed S = 1 Heisenberg chain is 
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bounded around the center of the system. In this section 
we focus on the excitation energy 



A(A) = 2^ (A) - E (X) 



(4.1) 



and investigate its dependence on A. Here, E (X) is the 
ground-state energy and E l (A) is the energy of the first 
excited state. As it is shown in the following, A(A) is 
insensitive to the boundary interaction parameter J cn d- 

4-1 Insensitivity of A(A) with respect to J en d 

Figure 8 shows A M (A, J cn d) f° r the 100-site system, 
where A M (A, J ond ) is the calculated gap for the M-site 
system when the boundary interaction is J ond . We have 
chosen the set of deformation parameter A = 0, 0.1, 0.2, 
0.3, and 0.4, with J ond = 0.25, 0.5, and 1. When A = 
0, the gap Aioo(0, Jend) is dependent on the value of 
J ond . This is because the excited quasiparticle can reach 
the system boundary, and it is affected by the effect of 
J ond . In particular, when A = and J cnd = 0.25 the 
quasiparticle is even localized near the system boundary. 
Appearance of this surface excitation tells that we need 
to avoid a parameter range of A <C 0.1 with J end = 0.25. 

On the other hand when A > 0.1, the excited quasi- 
particle cannot reach the system boundary as shown in 
Fig. 4, and the effect of J end on Ai o(A, J cn d) is negligible 
for this excitation. In this way, the hyperbolic deforma- 
tion has an effect of separating elementary excitations 
from the system boundary. 

To fix a desirable parameter range of M and A, in 
order to decouple quasiparticle from the boundary, let 
us observe the system size dependence of Am (A, J cn d) 
for the case A = 0.0025, which is the smallest one used 
in the following analysis. The system size is increased up 
to M = 2N = 1000. We introduce a quantity judging an 
error 

A M (A,J cnd ) 



£a = 



,(A,0.25) 



1 



(4.2) 



In the parameter region where e A is close to zero, we 
can say that the boundary effect is removed. Figure 9 
shows A M (0.0025, J ond ) with respect to 1/M = 1/2 AT 



at J = 0.25, 0.5 and 1. When J 



end 



0.25 and 



Fig. 9. Excitation energy Am(0.0025, J C nd) as a function of 
1/M = 1/(2JV) under the condition J end = 0.25, 0.5, and 1. 
The inset shows £a m eq. (4.2). 



when the system size is small, the boundary excita- 
tion is detected as appearance of smaller gap around 
A M (0.0025,0.25) ~ 0.26. Because the strength of the 
boundary interaction is J cnd coshAA, the energy of the 
boundary excitation increases with the system size. Fi- 
nally the bulk excitation A(A = 0.0025) is detected as a 
gap of Am (0.0025, Jcnd) — 0.41 in the neighborhood of 
M = 2N = 1000, From e A shown in the inset, we can 
say that the effect of J ond to Am (A, J cn d) is less than 
10~ 9 when the system size reaches 1000. Note that e A is 
greater than 10~ 10 since we set the convergence threshold 
of the eigensolver, which is used in the finite-size sweep- 
ing process, to the value 10~ 8 . When it is necessary, we 
decrease the threshold down to 10~ 10 in the following 
numerical calculations. 



4-2 A dependence of the energy gap 

We have erased the effect of system boundary from 
the elementary excitation by the hyperbolic deformation. 
Thus the estimation process for the Haldane gap A is re- 
duced to the extrapolation of A (A) with respect to the 
deformation parameter A. As we will see, the difference 
A(A) — A(0) is nearly proportional to A, where the de- 
pendence is consistent with the picture of the shallow 
bound state appearing in §3. 

In the following analysis, we use A 6 oo(A, 1) as the bulk 
excitation A(A), because the estimation value Am(A, 1) 
is not changed within numerical precision we require 
when the system size M runs from 600 to 1000. Fit- 
ting 2nd-order polynomials to the shown data, we ob- 
tain A(0) = 0.41047941(1) when m = 160 and A(0) = 
0.41047931(1) when m = 180, where the numbers in 
the parenthesis represent the fitting error. As we have 
discussed, we may choose the smaller one A(0) = 
0.41047931 as the estimated upper bound of the Haldane 
gap A. 

Let us check the precision of the value A(0) = 
0.41047931(1) that we have obtained from the indepen- 
dent data analysis. Let us assume that A(A) can be rep- 
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Fig. 11. A(A',<5) in eq. (4.5) with respect to A'. The dotted line 
represents the linear fitting. 



resented in terms of a polynomial in A 

A(A) = A(0) + aX + bX 2 + 0(A 3 ) . (4.3) 

A way of estimating A(0) efficiently is to consider the 
derivative between two values of deformation parameters 

A 2 A(A 1 )-A 1 A(A 2 ) 



X 2 -X 1 



that the method enables one to estimate a lower (upper) 
bound simultaneously with an upper (lower) bound of 
the gap. 

Here, we review this estimation method. First, we con- 
sider the initial sequence A$ for M = 2,4,6,..., which 

(fe) 

M 



— A(0) — b\ 1 X 2 + O (A 3 ) , (4.4) i s convergent to ■ Suppose that we can generate A 



which does not contain first order term in A. Introducing 
the notation A' = (A 2 + A x ) /2 and 5 = (A 2 — A x ) /2, this 
derivative can be written as 

A(0) - b (A' 2 - <S 2 ) + O (A' 3 ) = A (A', 6) . (4.5) 

Using the data shown in Fig. 10, we calculate A(A',<5) 
for neighboring As and plot the result in Fig. 11. It is 
obvious that the third order correction 0(A' 3 ) is very 
small in the shown parameter area, and A(X',S) is al- 
most linear in A' — S 2 . By use of linear extrapolation we 
obtain A(0) = 0.41047940(2) when m = 160 and A(0) = 
0.41047931(1) when m = 180, where we have shown the 
fitting error in the parenthesis. As we have done for the 
previous estimation for the upper bound, considering in- 
dependence of the value, we choose A(0) = 0.41047931 
as the candidate. This value is the same as the number 
obtained in the previous paragraph. 

5. Gap Estimation by the Sequence Interval 
Squeeze Method 

In this section, we present another new result on upper 
and lower bounds of the true Haldane gap using an ex- 
amination of numerical-diagonalization data. The reason 
for the usage of the numerical-diagonalization technique 
is to keep the high precision in the original numerical 
data. Having an independent estimate by a complemen- 
tary approach to the DMRG calculation, we obtain def- 
inite values on two bounds, which are concluded in the 
last section. 

Quite recently Nakano, one of the authors, and Terai 
proposed a new way to create an increasing (decreas- 
ing) sequence from monotonically decreasing (increas- 
ing) sequences of numerical data, which was immedi- 
ately applied to estimation of the Haldane gap with high 
accuracy. 8 A noticeable superiority of their analysis is 



for k = 1,2,3,... from Ay by some method so that 
the new sequences A M ^ have the same limit, namely 
A^ = A^ for all k. Supposed also that the set of the 
sequences A^ has the following properties, 

(i) A^j is monotonic with respect to M. 



(ii) £m increases with M ; > £ (fe) 

(iii) decreases with fc; £m +1 ^ < Cm ■ 
where £ is the decay length 7 given by 



M ! 

(fe) 



/ A (k) - A { 

Cm = 2 / log 



,(fc) 



,(fc) 

l M-2 



l(*0 



(5.1) 



Nakano and Terai introduced another sequence 8 ob- 
tained from the sequences A M ^ and A^ which are con- 
vergent to A&$ from the same side. The new sequence 
(fe) 

By is given by 



B 



A (k) Ak') _ Ak) .(k 1 ) 
M Ji M+2 I± M 



(k) 

M + 1 Ak') _ Ak') _ .(k) Ak) ' 



(5.2) 



(fe) 

for k > k . The most important property of B M is 
that the new sequence is convergent to the same limit 
A& from the opposite side. (See appendix in ref. [8].) 
Thus, there is a relation B M \ ] < A^ < A^J (or 



A 



\i < A^ < B^J) and we obtain a reliable in- 
terval including the limit A^ , which we would like to 
know. Here ki, k 2 , Mi and M 2 are integers. When the 
sequences A M and A M •* are monotonically decreasing, 
Txdnk,M(A]tf) and maxk,M(B^) are an upper bound 
and a lower one for A^ , respectively. When the direc- 
tion of Aj$ and A M ^ is opposite, maxj^Af (Aj^) and 
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min/5 ; M(-B^) give a lower bound and an upper one, re- 
spectively. 

(k) 

An appropriate set of the sequences A M ' can be gener- 
ated systematically from A$ to make the interval nar- 
rower by using convergence-acceleration techniques. As 
such a technique, we discuss the e-algorithm 7 ' 23 and its 
generalization. The e-algorithm provides us with a new 
sequence of one-level higher, by the relation between the 
neighboring three levels 

1 (5-3) 



Table II. Sequence of finite-size gaps by cq. (5.3) with a = 0.4. 



i(fc+l) 
1 M 



1 M-2 



1 



A (k) 



A (k) 
■"-M-2 



+ 



1 



a 



A {k) 

M 



A(k-l) 



A (k) ■ 



for a = 1. This algorithm for a = 1 was originally de- 
veloped by Wynn. 23 In order to create the first level se- 
quence A$ , we can prepare a dummy sequence A^ M ^ 
of the level —1, where all of the elements are infinitely 
large, namely A^ 1 ^ = oo. 

To know whether the acceleration is successful or not, 
the monitoring of the decay length given by eq. (5.1) 
is important. One can consider that the convergence of 
a transformed sequence A$ for k > 1 is successfully 
accelerated when all the above conditions (i)-(iii) hold. 

We put the excitation gap of the M-sitc S = 1 chain 
Tt = X)i=i • Si + i under a certain boundary condition 
into the initial sequence A$ . Note that the direction of 
the monotonic behavior of the sequence depends on the 
boundary condition adopted. If we impose the periodic 
(S M +i = Si) or the twisted (<Sm+i = ~Sf,S y M+l = 
-Sf , Sm+i = Sf) boundary condition, the direction is 
different with each other. 8 What we would like to do 
is to estimate A^ only from a finite part of the initial 
sequence. 

For long, only the systems under the periodic bound- 
ary condition were examined in most of the finite-size- 
scaling studies based on numerical diagonalization data. 
Finite-size gap of such systems usually decreases mono- 
tonically, when the system size is increased. If we ap- 
ply the e-algorithm to the monotonically decreasing se- 
quence, we obtain only upper bounds of the Haldane gap. 
Nakano and Terai found that the excitation gap of the 
finite size systems under the twisted boundary condition 
is monotonically increasing. 8 Thus, we should examine 
which is better among both of the boundary conditions, 
periodic and twisted. 

We can interpret that Nakano and Terai found a quite 
systematic approach to overcome limitation of the usage 
of the original e-algorithm. 8 Fundamental steps are sum- 
marized as usage of 1) examination of different boundary 
conditions, 2) acceleration of the monotonic sequences 

Ajj, and 3) the above new sequence B$ . If an appro- 
priate acceleration transformation is chosen, one easily 
obtains a reliable interval including the limit Ao$ , which 
we would like to know. The interval gets narrower as 
the number of initial data is increased. Hereafter, we call 
the above procedure the sequence interval squeeze (SIS) 
method. Getting the bounds of both sides, we can quan- 
titatively discuss the precision of the estimates of the 
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4 (2) 
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AV 


A$ 

M 


AS) 
>M 


12 


0.409599020 








14 


0.410158700 








16 


0.410354714 


1.91 


0.410442366 




18 


0.410427448 


2.02 


0.410465180 




20 


0.410456146 


2.15 


0.410473181 


1.91 


22 


0.410468261 


2.32 


0.410476497 


2.27 


24 


0.410473733 


2.52 


0.410477982 


2.49 


M 






4(5) 
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20 


0.410476976 








22 


0.410478577 








24 


0.410479051 


1.64 


0.410479218 





Table III. Antimonotonic sequence B^> for a = 0.4 



M 


o 2 

M 






13 


0.411197171 






15 


0.410724641 






17 


0.410568330 


0.410491523 




19 


0.410513257 


0.410482576 




21 


0.410492985 


0.410480499 


0.410480148 


23 


0.410485151 


0.410479829 


0.410479554 



Haldane gap that were reported so far. 

Nakano and Terai substituted the excitation gap of the 
finite size systems under the twisted boundary condition 
for A^j ; the initial sequence is monotonically increasing. 



From the numerical data collected up to A^ =2i with the 
twisted boundary condition, the SIS method provides us 
with the estimate A = 0.4104789(13) as the Haldane gap 
from the acceleration process under a = 1. 

Let us recall that an available acceleration transfor- 
mation is not limited to the e-algorithm with a = 1. 
Actually, the first step of e-algorithm is equivalent to 
the Aitken-Shanks transformation, 24 which corresponds 
to a = 0. If one takes a = for the second step and the 
later, the transformation (5.3) is reduced to just an itera- 
tion of the Aitken-Shanks transformation. Other choices 
of a for the second step and the later correspond to differ- 
ent convergence acceleration transformations, on each of 
which the degree of acceleration depends. Thus we may 
optimize the convergence acceleration in the step 2) in 
the SIS method. If we adjust a within a range in which 
the above three conditions are certified, the strength of 
the acceleration can be optimized. In the present work, 
we apply a = 0.4 in cq. (5.3) giving A$ and k' = k — 2 
in eq. (5.2) giving to the same numerical data with 
the twisted boundary condition. The result successfully 
gives an inside interval narrower than the result obtained 
by a = 1. 

Let us list all the A$ created by eq. (5.3) in Ta- 
ble II. Sequences A^j and A$ are not presented be- 
cause they do not change due to a variance of a. As it 
is observed, A$ and £^ satisfy the set of the three 
conditions for the successful convergence acceleration. 
We then create another sequence by eq. (5.2) and 
show them in Table III. One can clearly observe that 
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B ( M ' is monotonically convergent from the side oppo- 
site to Afy[ . Figure 12 depicts the 1/M dependence of 
these sequences. The present result gives a new interval 
[0.410479051, 0.410479554] for an estimate for the Hal- 
dane gap A, which is narrower than the reported one in 
ref. [8]. 

6. Conclusions and Discussions 

We have observed the magnetic excitation of the 
5 = 1 Heisenberg chain, whose exchange coupling is 
deformed hyperbolically. The magnetic quasiparticle is 
weakly bounded in the neighborhood of the center of the 
system. It is shown that the excitation energy A(A) is 
nearly linear in A, and the extrapolation to A = gives 
an estimate of the Haldane gap, as precise as that ob- 
tained by the SIS method. Since the quasiparticle does 
not reach the system boundary under the hyperbolic de- 
formation, one does not have to pay special attention to 
the boundary interaction strength. 

Here, let us summarize the estimated values of the 
Haldane gap we have obtained with those reported so 
far. Figure 13 shows A(0) obtained by 

A. Lanczos method + convergence acceleration in 
ref. [25] 

B. DMRG applied to undeformed system in ref. [13] 

C. Monte Carlo simulation in ref. [26] 

D. Lanczos method + SIS in ref. [8] 

E. DMRG applied to undeformed system in §3. 

F. Deformation analysis in §4 with polynomial fitting. 

G. Deformation analysis in §4 with derivatives. 

H. Lanczos method + SIS in §5. 

Each method among E at J on( j = 0.5088, F, G has four 
data points when m =120, 140, 160 and 180. In addition, 
there are two data points under the conditions m = 160 
and 180 for the method E at J cnd = 0.50866. 

In Fig. 13, A(0) given by one of methods E-G always 
monotonically decreases, when m is increased. We know 



Fig. 13. The values of the Haldane Gap which have been reported 
and that calculated in this article. A scale of the perpendicular 
axis is adjusted so that all error bar are shown in the inset. 

that the energy lift by the cut-off effect in the excited 
state is larger than that in the ground state. Besides, 
the systems treated in §4 have enough sizes, where any 
boundary effects are eliminated. Thus, we may suppose 
that the estimated values of A(0) converge monotoni- 
cally with respect to m > 120. Note that the differ- 
ence between A(0) = 0.41047941(1) when m = 160 and 
A(0) = 0.41047931(1) when m = 180 is big enough com- 
pared with the fitting error. 

It should be noted that the SIS estimation 8 gives both 
the lower and the upper bound for the Haldane gap. The 
values given by the item H provide us the best bounds 
with this method. When we choose m = 180, our esti- 
mations in E-G are always inside the interval given by 
the item H. 

Thus, we may safely conclude that 0.41047931 should 
be a better upper bound of the Haldane gap than 
the value given by the item H. Following this dis- 
cussion, we conclude that the Haldane gap is in 
[0.41047905, 0.41047931]. Looking at Figure 13, we may 
also suppose that the actual Haldane gap is closer to the 
upper bound than the lower bound. 

For the hyperbolic deformation, the choice of m is the 
remaining single parameter determining the accuracy. If 
we consider the m dependence of the gap, we may also 
construct a sequence interval squeeze technique by gener- 
alizing the method given in §5. We can thus find that the 
use of the hyperbolic deformation is one of the efficient 
tool to detect the excitation gaps of one-dimensional 
quantum systems. 

One of future subjects is numerical gap estimation of 
5 = 2 antiferromagnetic Heisenberg chain. Let us con- 
sider appropriate boundary conditions based on the va- 
lence bond solid (VBS) picture for this case. A possible 
simple choice is to put 5 — 1 spins at the both ends of 
the system. Another simple choice is to put two 5 = 1/2 
spins at the each end of the system, where the bond 
configuration is the form of the letter Y. In addition, a 
slightly complex choice is to reduce the length of spin by 
amount of 1/2 site by site, i.e. to put 5* = 3/2, 5 = 1, 
and 5 = 1/2 spins at the boundary. For each candidate of 
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the boundary spin arrangement, one has to consider the 
parametrization of the bond strength. In this way, the 
tuning of the boundary condition for the 5 = 2 chain is 
more complicated than that of S = 1 chain. Therefore, 
an efficient arrangement of additional boundary spins has 
not been reported. When the system is deformed hyper- 
bolically, however, the problem of parameterization could 
be put in a extrapolation of A. 

Another one is to find out classical analogue of the hy- 
perbolic deformation for 2D statistical models. A candi- 
date is the hyperbolic lattice models studied so far, 27-30 
but in those models only discrete values of A are allowed. 
To construct a class of models that has appropriate struc- 
ture for the DMRG applied to classical systems 31 would 
be important for the further study. 
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